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Abstract 

We present an algorithm which allows a fast numerical computation of Feldman-Cousins 
confidence intervals for Poisson processes, even when the number of background events is 
relatively large. This algorithm incorporates an appropriate treatment of the singularities 
that arise as a consequence of the discreteness of the variable. 
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1 Introduction 

In physics there are many situations where the outcome of an experiment is a positive integer 
number with a Poisson distribution. This is the case for instance of the number of events of a 
certain type produced in high energy collisions. The statistical analysis of these processes is a 
difficult task when the result obtained is in the limit of the sensitivity of the experiment. In 
general, the number of events no obtained in an experiment consists of background events with 
known mean b and signal events, whose mean \x is the quantity that we want to determine. The 
problem arises when the number of events obtained no is significantly lower than the background 
expected b. This happens in some experiments on neutrino oscillations, for instance in the 
KARMEN 2 experiment 0. 

Usually, after performing an experiment, one decides whether to give the results on the un- 
known parameter \x in the form of a central confidence interval or an upper bound. This decision 
(called 'flip-flopping') is based on the data and, as has been shown by Feldman and Cousins [g], 
introduces a bias that may cause that the intervals cover the true value [i with a smaller frequency 
than the stated confidence level. To solve this and other problems, they introduce a new ordering 
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principle that unifies the treatment of central confidence intervals and upper limits. This is possi- 
ble because the Neyman construction of confidence intervals || allows the choice of the ordering 
principle with which the intervals are constructed. Typical choices lead to the construction of 
either central intervals or upper confidence limits. The choice of Feldman and Cousins gives inter- 
vals that are two-sided or upper limits depending on the result of the experiment and not on the 
choice of the experimentalist. These intervals avoid the undercoverage caused by 'flip-flopping' 
and are non-empty in all cases. Some variants of their method have been also proposed |4], B. 

To consider the Feldman-Cousins confidence intervals as an alternative to standard intervals, 
in practice one needs to calculate these intervals for arbitrary uq and b. The Tables provided 
in Ref. 0, for no < 20, b < 15, are sufficient for small luminosity /statistics experiments, but 
for higher luminosities in general b and no are larger. One possibility is to extend these Tables 
using the same systematic computational method of Ref. whose speed is not optimized and 
consumes a lot of time. More convenient is to develop a program which takes no and b as inputs 
and gives as output the confidence interval, requiring a minimal number of calculations. This is 
what is done here. The program can be used either directly to compute the confidence interval 
for given no and b or in conjunction with other routines. This is especially useful, for instance 
to calculate expected limits from rare high energy processes for different values of the center of 
mass energy or the collider luminosity, which is the case that we were primarily interested in. 

In the following we introduce a procedure to compute the Feldman-Cousins intervals in an 
efficient way for arbitrary no and b, in principle only limited by the machine precision. In Section 
2 we review Neyman's construction of the confidence intervals for a Poisson variable, emphasizing 
some points that simplify the numerical calculation. Section 3 is more technical and devoted to 
explain in depth how to translate this method for the computer calculation. In Section 4 we 
present our results. The FORTRAN implementation of the algorithm is given in the Appendix. 
Other implementations in C and Mathematica Q (about 100 times slower than the FORTRAN 
version) can be obtained from the author. 

2 Construction of the confidence intervals 

The probability to observe n events in a Poisson process consisting of signal events with unknown 
mean [i and background events with known mean b is given by the formula 

P(n|u;6) = ^±^e-^ 6 ), (1) 

with n, fj,,b > 0, and n restricted to integer values. The construction of the confidence intervals 
on the unknown variable fi follows Neyman's method of the confidence belts. 

The first step in this procedure is to construct, for a fixed value of b and for different values of 
/i, the confidence intervals [ni(/i; b),ri2(n; b)] such that the probability to obtain a result between 
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n\ and n 2 is greater or equal than a, the confidence level (C. L.), 
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P(ne \n u n 2 )\n;b) = £ P(n | /x; b) > a . (2) 

n=ni 

It is worth to note that for the more general case of a continuous variable x the intervals 
[xi(/j,; b), a^O-t; b)] satisfy P{n e [2:1,2:2] I A 4 ! b) = a. For a discrete variable n it is not possible 
to obtain the exact equality, and to avoid undercoverage it is replaced by the inequality in Eq. 

is. 

The choice of the intervals [ni,n2] is not unique, and determines the type of confidence 
intervals on /i that are constructed. The most common choices are n 2 = 00, P(n < n\ \ fx; b) < 
1 — a, which gives upper confidence bounds, and P(n < <n\ \ /x; b) < (1 — a)/2, P(n > n 2 \ /J,;b) < 
(1 — a)/2 which leads to central confidence intervals. The prescription of Ref. is based on a 
likelihood ratio R, constructed as follows. 

1. For any values of b and n, one considers which value of [i would maximize the probability 
P(n I fi; b) . It is straightforward to find that for n > 1, P(n \ fx; b) considered as a function 
of jj, grows for [i < n — b, has a maximum at /i = n — b and decreases for /i > n — b. As /i is 
restricted to lie in the positive real axis, if n > b the maximum is fi = n — b, otherwise the 
maximum is /i = 0. In the case n = the maximum is also /z = 0, so we define 

Mbest(«; b) = max{0, n-b} (3) 

as the value which maximizes P{n \ /x; b). 

2. Then, for any value of /i we consider the quantity R(n; b) defined as 

„/ , v Pin \ix\V) 

R(n;n,b)= y 7 J <1, 4 

on which the Feldman-Cousins ordering principle is based. To construct the interval 
[wi(/x; b), ri2(fJ,; b)], for each value of [i (and fixed b) one takes values of n with decreas- 
ing R(n;n,b), summing up their probabilities P(n\fi;b) until the total equals or exceeds 
the C. L. desired a. Thus the interval [ni,ri2] is the set of values of n necessary to satisfy 
the inequality in Eq. @, taken with the largest R(n; /1, b). 

The simplicity of this prescription allows a fast computer implementation. Instead of generating 
a large table of values for n and taking those with the largest R, we can directly find these values 
and add them successively to construct the interval. (This is difficult to do with the more involved 
prescriptions of Refs. [Q, ||.) For this purpose we will examine the behaviour of R. If we consider 

_|_ fyx e -(v+b) 

R{X] ^ b) = (f, hest (x;b)+bre-(^^b)+b) ( 5 ) 
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as a function of the continous variable x, we can look for its maximum. Let us first consider 
\i 7^ 0, 6 ^ 0. For x sufficiently small, Hhest(%] 6) = and R(x; /i, b) = (1 + /i/6) x e _/i is increasing. 
For x > b, fibest(x; b) = x — b and R(x; fi, b) = (e(/U + b) / x)) x e~^ +b ^ grows for x < fi + b, falls for 
x > fi + b and has a local maximum at x = fi + b, which is then the global maximum. This is also 
true for \i ^ 0, b = 0. For fj, = 0, b ^ and x < b, n^ es t(x; b) = and i?(x; 0, b) = 1, its maximum 
possible value. For x > b, R(x;0,b) = (b/x) x e x ~ h decreases with x. Hence the maximum is still 
x = fj, + b, although not unique. The only remaining case with \i = 0, b = in which the Poisson 
distribution is singular must be treated separately. 
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Figure 1: Confidence belt for 6 = 3, with a = 0.9. The vertical line is drawn for no = 8. 

With this method, and for different values of (b is fixed), one calculates the confidence 
intervals [m(/z; 6), ri2(//; £>)] obtaining a confidence belt like the one showed in Fig. |l] for 6 = 3 
and a C. L. of 0.9. To find the confidence interval [/ii,//2] for a particular experimental value no, 
one draws a vertical line at n = no and finds the maximum and minimum values of \i for which 
the line intersects the confidence belt. In Fig. |l| we observe that [i\ is the smallest /i such that 
"2(a*; b) = no, whereas fj>2 is the largest /i for which ni(/x; 6) = no- 

3 Implementation 

Let us explain how the method described in Section 2 is made suitable for the evaluation in a 
computer. For the calculations we use the FORTRAN version of the program compiled with f ort77 
under Linux on a Pentium III-450 (compiled with g77 the program runs about 25% slower), and 
for the plots we also use the Mathematica version. 




4 



The first problem in the practical realization of the Neyman construction is that, for large 
b, n or /x, the factors of the Poisson probability formula in Eq. ([l]) can overflow (or underflow) 
the computer capacity in intermediate calculations. The factor (/i + b) n may be very large, for 
instance 54 178 is larger the biggest double precision real number in the FORTRAN compiler used, 
approximately 10 308 . However, the exponential factor in Eq. ([l]) compensates for it in the final 
result. For \x + b + n > 230 we evaluate P using the expression 

P(n\^b) = e-»f[(^)e- b , (6) 

i=i v 1 J 

with the product calculated factor by factor. This extends the allowed size of the parameters 
of our program, with the disadvantage of a larger computing time. For \i + b + n < 230 the 
factor (/j, + b) n is not too large, and can be directly calculated. In this case the factorials up to 
170! ~ 10 are calculated at the beginning of the main program and stored in the array fact to 
save time, whereas for n > 170 the expression of the Poisson formula is divided by n£=i i factor 
by factor. 

The quantity R in Eq. @ is a ratio of probabilities and can be computed without any problem 
cancelling out the common factors and defining a function R. (Defining /ibest as DIM (FLOAT (n) ,b) 
instead of MAX (OdO, FLOAT (n)-b) as we do would not have improved the speed significantly.) 

The core of the algorithm is the subroutine NRANGE, used to calculate the confidence intervals 
[ni(fj,; b), ri2 (a*; b)] for arbitrary \i and b. Its arguments are the variables rmu (//), b and the 
confidence level desired CL. The output nl and n2 is given in a COMMON block, together with a 
variable CLac, the C. L. finally achieved (in general it is greater than CL) which is useful for 
other purposes. The discussion in the last Section simplifies the implementation of the algorithm 
considerably, because we have found that the values of n that maximize R(n; /i, b) concentrate 
around fi+b. This improves the speed by an order of magnitude for large values of the parameters, 
since we do not need to calculate a large table (n, P(n), R(n)) and sort it. Instead, we know the 
maximum R is one of the two integers nearest to /i + b. We begin with nl=INT(rmu+b) , n2=nl+l. 
If R(nl,rmu,b) is larger than R(n2,rmu,b) we take nl, decrease nl and add P(nl,rmu,b) to 
CLac. Otherwise, we take n2, increase n2 and add P(n2,rmu,b) to CLac. Repeating this until 
CLac is greater than CL and taking into account that nl must be greater than zero we obtain the 
desired interval [nx, 77.2]- The singular case /i = b = is treated separately. With the Mathematica 
version of this subroutine we can plot confidence belts like that in Fig. |l[ 

The calculation of the confidence interval [^1,^X2] is done using two funcions RMUl(n,b,CL) 
and RMU2(n,b,CL), where n is the experimental number of events. We discuss them in turn. 

As we see in Fig. |], the lower limit [i\ is the minimum value of /i such that 77-2 (yu; 6) = 
riQ. Within our framework, the calculation is done looking for the minimum rmu such that n2 
calculated with NRANGE (rmu, b,CL) equals n. The search is done with the bisection method. 
Starting with the limits rmumin=0d0, rmumax=FL0AT(n) -b+ldO (RMU1 must be between these two 
values) we calculate the midpoint of the interval, rmumed, and NRANGE (rmumed,b,CL). If n2 is 
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greater or equal than n, we move rmumax to rmumed, otherwise we move rmumin to rmumed. This 
is repeated until the length of the interval is smaller than the desired precision delta, which we 
take as the maximum of 0.01 and 0.0005 times the background b. We summarize the algorithm 
in Fig. g 

Set initial 
mumin and mumax 





New mumax 



Figure 2: Flux diagram for the calculation of RMU1. 



In principle the calculation of RMU2 would follow an analogous procedure. However, in this 
case we find an extra problem. Except for a few cases with small 6, the confidence belt is not as 
simple as in Fig. [l] but displays a more elaborated structure as can be seen in the example of 
Fig. |. 

The fact that ni(u;6) is not a monotonic function of u is due to the discreteness of n. This 
causes that the set of u values for which the vertical line at n = no intersects the belt is not 
connected for no = — 4. The effect is relevant since the upper limit U2 is defined as the largest 
\i for which ni(u; b) = uq. Thus a modification of the algorithm is required not to miss the small 
wedges in the function. 

For n > b the behaviour is as expected and we can use the same algorithm as for fj,\. We 
have checked values of b between and 50 and have found that for n > b the function ni (//;&) 
does not have any singularity, so in this case we can safely adapt the routine RMU1. We look for 
the maximum rmu such that nl calculated with NRANGE(rmu,b,CL) equals n. The search is again 
done with the bisection method. We start with the initial values rmumin=MAX (OdO, FLOAT (n)-b), 
rmumax=3dO*SQRT (FLOAT (n)+b+ld0) . (If rmumax is not sufficiently large, we increment it in steps 
of SQRT (FLOAT (n) +b+ld0) .) We calculate the midpoint of the interval, rmumed, and calculate 
NRANGE(rmumed,b,CL) . If nl is lower or equal to n, we move rmumin to rmumed, otherwise we 
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Figure 3: Confidence belt for b = 7, with a = 0.95. 

move rmumax to rmumed. This is repeated until the length of the interval is smaller than the 
desired presision delta. 

For n < b we sample the interval for possible singularities, which consumes more time. 
This is done in three iterations with increasing number of points. To minimize the length of 
the interval and optimize the density of the sampling, we take rmumin=MAX(OdO, FLOAT (n)-b) 
and increase it in steps of 1 while nl is lower or equal to n. As the upper limit we take 
rmumax=3dO*SQRT (FLOAT (n) +b+ldO) , sufficiently high so that the initial interval contains all 
the singularities for a C. L. of 0.99 or less. 

In the first step we divide the interval between rmumin and rmumax in 10 parts and check if 
any of the points selected has nl lower or equal to n. If it is so, we change rmumin to the largest 
of them (this is always safe) and start again with this new rmumin. This first sampling with a 
small number of points finds wedges like those for n = 1, 3 in Fig. ^ and saves a lot of computing 
time. The narrow wedges at n = 0, 2 require more dense samplings. 

If the points calculated have nl greater than n, we check the upper half of the interval for 
singularities. The second iteration divides the upper half in 20 parts and checks 19 points. If it 
finds any singular point with nl lower or equal to n, it changes rmumin and starts again at the 
first step. If not, the third iteration divides the upper half in 500 parts. If a singular point is 
found, it changes rmumin and starts the first step. If not, unless some kind of singular behavior 
is found (in which case a fourth sampling with 5000 points is performed) it is assumed that there 
do not exist singular iries and rmumax is changed. 

An additional speed improvement is implemented: if in the second or third iterations the 
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density of points is sufficiently high (the points are closer than stepmin), the number of points 
is decreased and no more iterations are performed. The flux diagram of RMU2 is shown in Fig. |j. 
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Figure 4: Flux diagram for the calculation of RMU2. 



To simplify changing the parameters of this routine, the number of iterations and their re- 
spective number of points are stored in the variables maxit and maxdivs. The calculation for 
n > b is done with the same function with maxdivs (1) =2 and only the first iteration. 

One may notice that some upper limits obtained with RMU2 are different from those quoted 
in Ref. Qj. This is again a consequence of the discreteness of n. The upper limit \xi for no fixed 
is not always a decreasing function of b (dotted lines in Fig. ||.) This behaviour is corrected in 
Ref. ||] forcing the function to be nonincreasing, calculating the upper limit \i2 from b = 25 to 
b = in steps of —0.001. The corrected value is then the maximum of fj,2(no,x) for x > b. Some 
people, however, find this ad hoc correction questionable Q. At any rate we could also follow 
this procedure getting the same values of Ref. and the solid line in Fig. ||. This requires a 
very long calculation (25000 different values of b), for instance the time to calculate the uq = 
line is 34 m 27 s. 

Of course, to obtain \ii for a particular b it is not necessary to calculate the whole interval 
[b, 20] and it is enough to consider approximately [b, 6+1]. For this purpose we use the function 
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RMU2c. This routine examines the behaviour of RMU2 in the interval [b,b + 1] and corrects the 
value if necessary. The adjacent maxima that can be seen in Fig. [| are found with the simple 
golden section search of Ref. |7]]. (Other more sophisticated methods offer no advantage since 
the function does not seem to be differentiable at the maximum.) The initial bracketing of the 
maximum is very delicate as can be also seen in this Figure. We do it examining RMU2(n,bl,CL) 
taking bl with increments of 0.1 until RMU2 begins to grow, then in increments of 0.05 until 
it begins to decrease. Then the golden section method is applied to find the maximum with a 
precision of 0.001. This maximum is then compared to the value at b to take the largest value. 

This method again brings a substantial speed improvement over the blind computation in 
steps of —0.001 in 6, as we will see in next Section, Examining the behaviour of (12 we have also 
found that the correction is not necessary in general for n > b/2 and the function RMU2 could be 
used directly. There are however some exceptions, for instance n = 10, b = 14 with a C. L. of 
0.95. To be conservative, we will only use RMU2 when n > b. 

4 Results 

To obtain our results we use the same precision that is used in Ref. ||, a minimum step stepmin 
of 0.005 in /i and an accuracy of 0.01 in the upper and lower limits of the confidence intervals. 
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To calculate the limits in their Tables II-IX including the singular cases it is enough to con- 
sider in the third iteration maxdivs(3)=100 for confidence levels of 68.7%, 90% and 95%, and 
maxdivs(3)=300 for a C. L. of 99%. For better comparison we use maxdivs(3)=500 as we do 
in the rest of the calculations to ensure that all singularities are found. The running time is 
summarized in Table [l|. 

C. L. h t 2 
68.7% 4.1 s 1 m 9.0 s 
90% 6.1 s 1 m 28.6 s 
95% 6.8 s 1 m 36.8 s 
7.9 s 1 m 57.4 s 



Table 1: Time spent in the calculation of the confidence intervals for tiq < 20 and b < 15, with 
RMU2 (ti) and with RMU2c (t 2 ). 



To check if our algorithm in fact handles large numbers efficiently we measure the time spent 
to calculate the upper limit [i 2 with RMU2 for hq = b between and 200, obtaining the solid line 
in Fig. |(| We can also use RMU2c forcing the program to look for unexisting spurious maxima in 
b (and hence also for singularities with maxdivs(3)=500) obtaining the dotted line. 
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Figure 6: Time required to calculate the upper limit fi 2 with the function RMU2 (ti) and with 
RMU2c forcing to look for adjacent maxima (t 2 ) as explained in the text. 
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It is amazing to observe that the computing time not only does not grow quickly with b as 
it could be expected, but remains almost constant for b < 75. This is achieved with (i) a fast 
algorithm to find the singularities if they exist, (u) the optimization of NRANGE to calculate only 
the data really needed, and (Hi) the calculation of the factorials up to 170! at the beginning of 
the program. For b > 75 rmu+b+FLOAT (n) is sometimes larger than 230 and the time required 
begins to grow linearly with b after the gap between 75 and 100, as can be observed in the Figure. 
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A Appendix 

PROGRAM DEMO 

IMPLICIT REAL* 8 (A-H.O-Z) 
DIMENSION FACT (0:170) 
COMMON /range/ nl,n2,CLac 
COMMON /factorial/ FACT 
DATA CL /0.99d0/ 

FACT(0)=ld0 ! FACT initialization 

DO i=l,170 

FACT ( i ) =FACT ( i - 1 ) *DFL0 AT ( i ) 
ENDD0 

C CALCULATION OF TABLES VIII, IX OF REF [2] 

DO b=0d0,4d0,0.5d0 
DO n=0,20 

PRINT 100,n,b,RMUl(n,b,CL) ,RMU2c(n,b,CL) 
ENDD0 
ENDD0 

DO b=5d0,15d0,ld0 
DO n=0,20 

PRINT lOO.n.b.RMUKn.b.CL) ,RMU2c(n,b,CL) 
ENDD0 
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ENDDD 
STOP 

100 FORMAT (>n = ',12,', b = ',F4.1,' -> [ ' ,F6 . 2 , ' , ' ,F6 . 2 , ' ] ' ) 
END 



DOUBLE PRECISION FUNCTION P(n,rmu,b) 
IMPLICIT REAL* 8 (A-H.O-Z) 
DIMENSION FACT (0:170) 
COMMON /factorial/ FACT 

IF ((rmu .EQ. OdO) .AND. (b .EQ. OdO)) THEN 
P=0 

IF (n .EQ. 0) P=l 
RETURN 
END IF 

IF (rmu+b+FLOAT(n) .LE. 230d0) THEN 
P= (rmu+b) **n*EXP (-rmu-b) 
IF (n .LE. 170) THEN 

P=P/FACT(n) 
ELSE 

P=P/FACT(170) 

DO i=171, n 
P=P/FLOAT(i) 

ENDDO 
END IF 
ELSE 

P=EXP(-rmu) 
DO i=l,n 

P=P* (rmu+b) /FLOAT (i) 
ENDDO 

P=P*EXP(-b) 
END IF 
RETURN 
END 



These lines 
are not 

needed if P is 
only called 
from NRANGE 



This is not normally used because for 
n > 170 rmu+b+FLOAT(n) will be larger 
than 230 



DOUBLE PRECISION FUNCTION R(n,rmu,b) 
IMPLICIT REAL* 8 (A-H.O-Z) 
IF (n .LT. 0) THEN 
R=OdO 
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RETURN 
END IF 

R=EXP (MAX (OdO , FLOAT (n) -b) -mm) 

IF (n .GT. 0) R=R* ( (rmu+b) / (MAX (OdO , FLOAT (n) -b)+b) ) **n 

RETURN 

END 



SUBROUTINE NRANGE (mm , b , CL) 
IMPLICIT REAL* 8 (A-H.O-Z) 

COMMON /range/ nl,n2,CLac ! CLac for future use 
IF ((rmu .EQ. OdO) .AND. (b . EQ. OdO)) THEN 

nl=0d0 ! Special case 

n2=0d0 

RETURN 
END IF 

nl=INT (rmu+b) ! The maximum R is between 

n2=nl+l ! these values 

rl=R(nl ,rmu,b) 

r2=R(n2,rmu,b) 

CLac=0d0 

DO WHILE ((CLac .LT. CL) .AND. (nl .GE. 0)) 
IF (rl .GT. r2) THEN 

CLac=CLac+P (nl , rmu , b) 

nl=nl-l 

rl=R(nl ,rmu,b) 
ELSE 

CLac=CLac+P (n2 , rmu , b) 
n2=n2+l 
r2=R(n2,rmu,b) 
END IF 
ENDDO 

DO WHILE (CLac .LT. CL) ! No need to calculate R 

CLac=CLac+P (n2 , rmu , b) 

n2=n2+l 
ENDDO 
nl=nl+l 
n2=n2-l 
RETURN 
END 



13 



DOUBLE PRECISION FUNCTION RMUl(n,b,CL) 
IMPLICIT REAL* 8 (A-H.O-Z) 
COMMON /range/ nl,n2,CLac 
CALL NRANGE ( OdO , b , CL ) 
IF (n2 .GE. n) THEN 

RMUl=OdO ! Special case 

RETURN 
END IF 

rmumin=OdO 

rmumax=FLOAT(n) -b+ldO 
delta=MAX (0 . OldO , . 0005d0*b) 

DO WHILE ( (rmumax-rmumin) .GE. delta) ! Bisection method 

rmumed= (rmumin+rmumax) /2d0 
CALL NRANGE (rmumed,b,CL) 
IF (n2 .GE. n) THEN 

rmumax=rmumed 
ELSE 

rmumin=rmumed 
END IF 
ENDDO 

RMU1= (rmumin+rmumax) /2d0 

RETURN 

END 



DOUBLE PRECISION FUNCTION RMU2(n,b,CL) 

IMPLICIT REAL* 8 (A-H.O-Z) 

DIMENSION maxdivs(4) 

LOGICAL safe , saf enow, sing, changemin 

COMMON /range/ nl,n2,CLac 

DATA maxdivs /10 , 20 , 500 , 5000/ 

maxit=4 

stepmin=0.005d0 
rmumin=MAX (OdO , FLOAT (n) -b) 
CALL NRANGE ( rmumin , b , CL ) 
IF (FLOAT (n) .LT. b) THEN 

DO WHILE (nl .LE. n) ! Take lower limit 

rmumin=rmumin+ldO ! as high as possible 
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CALL NRANGE (rmumin , b , CL) 
ENDDD 

rmumin=rmumin- IdO 
safe=. FALSE. 
ELSE 

maxdivs(l)=2 ! Use bisection method when 

safe=.TRUE. ! there aren't sing. 

END IF 

rmumax=3dO*SQRT (FLOAT (n)+b+ldO) ! Large enough for most purposes 
CALL NRANGE ( rmumax , b , CL ) 

DO WHILE (nl .LE. n) ! If not, increase it 

rmumax=rmumax+SQRT(FLOAT(n)+b+ldO) 

CALL NRANGE (rmumax, b,CL) 
ENDDO 

delta=MAX (0 . OldO , . 0005d0*b) 
DO WHILE ( (rmumax-rmumin) .GE. delta) 
step= (rmumax-rmumin) /FLOAT (maxdivs ( 1 ) ) 
rmumin2=rmumin 
DO i=l,maxdivs(l)-l 

CALL NRANGE (rmumin+FLOAT ( i) *step , b , CL) 
IF (nl .LE. n) rmumin2=rmumin+FL0AT(i) *step 
ENDDO 

IF (rmumin2 .GT. rmumin) THEN 

rmumin=rmumin2 ! New rmumin -> change it 

ELSE 

safenow=safe ! Have to look for singularities 

sing=. FALSE. ! if they may exist 

changemin= . FALSE . 
it=2 

DO WHILE ((safenow . EQ. .FALSE.) .AND. (it .LE. maxit)) 
ndi vs=maxdi vs ( it ) 

step= (rmumax-rmumin) /FLOAT (2*ndivs) 

IF (step .LT. stepmin) THEN ! step is small 

ndivs=INT( (rmumax-rmumin) /(2*stepmin) )+l ! enough and this 
step=(rmumax-rmumin)/FL0AT(2*ndivs) ! will be the 

saf enow=.TRUE. ! last iteration 

END IF 

CALL NRANGE ( (rmumin+ rmumax) /2d0 , b , CL) 

n_prev=nl 

DO i=l,ndivs-l 
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CALL NRANGE (rmumin+FLOAT ( i+ndivs) *step , b , CL) 
IF (nl .LE. n) THEN 

rmumin2=rmumin+FLDAT( i+ndivs) *step ! New rmumin 

changemin= . TRUE . 
END IF 

IF (nl .LT. n_prev) sing=.TRUE. 
n_prev=nl 
ENDDO 

IF (changemin . EQ. .TRUE.) THEN 

rmumin=rmumin2 ! Change rmumin 

saf enow=.TRUE. ! Exit loop 

ELSE 

IF ((sing . EQ. .FALSE.) .AND. (it . EQ. maxit-1)) THEN 
safenow=.TRUE. ! Enough iterations 

END IF 
END IF 

it=it+l ! Next iteration 

ENDDO 

IF (changemin .EQ. .FALSE.) rmumax=(rmumin+rmumax)/2d0 
END IF 
ENDDO 

RMU2= (rmumin+rmumax) /2d0 

RETURN 

END 



DOUBLE PRECISION FUNCTION RMU2 
IMPLICIT REAL* 8 (A-H.O-Z) 
LOGICAL sing 

PARAMETER (R=0 . 61803399d0 , C=ld 

RMU2c=RMU2(n,b,CL) 

IF (n .GE. INT(b) ) RETURN 

stepl=0. IdO 

step2=0.05d0 

deltab=0.001d0 

blmax=b+ldO 

sing=. FALSE. 

C Bracket the maximum, if any 



(n.b.CL) 
-R) 

! Do not need correction 

! Go downhill in steps of 0.1 

! Go uphill in steps Of 0.05 

! Final precision in b 

! Look for maximum up to b+1 
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bl=b 
al=RMU2c 

DO WHILE ((bl .LE. blmax) .AND. 
a_next=RMU2(n,bl+stepl,CL) 
IF (a_next .GT. al) THEN 

sing=.TRUE. 
ELSE 

al=a_next 
bl=bl+stepl 
END IF 
ENDDD 

IF (sing .EQ. .FALSE.) RETURN 

b2=bl+stepl-step2 

a2=a_next 

DO WHILE (a_next .GE. a2) 
a2=a_next 
b2=b2+step2 

a_next=RMU2 (n , b2+st ep2 , CL) 
ENDDD 

b4=b2+step2 
a4=a_next 

IF (RMU2c .GT. a2+0.05d0) RETURN 
IF (b4-b2 .GT. b2-bl) THEN 

b3=b2+C*(b4-b2) 

a3=RMU2(n,b3,CL) 
ELSE 

b3=b2 

a3=a2 

b2=b3-C*(b3-bl) 
a2=RMU2(n,b2,CL) 
END IF 

C Find the maximum 

DO WHILE (b4-bl .GE. deltab) 
IF (a3 .GT. a2) THEN 
bl=b2 
b2=b3 

b3=R*b2+C*b4 
al=a2 



(sing .EQ . .FALSE.)) ! Go downhill 



! RMU2 is always decreasing 



! Go uphill 



! This maximum will not be larger 
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a2=a3 

a3=RMU2(n,b3,CL) 
ELSE 
b4=b3 
b3=b2 

b2=R*b3+C*bl 

a4=a3 

a3=a2 

a2=RMU2(n,b2,CL) 
END IF 
ENDDD 

RMU2c=MAX (RMU2c , a2 , a3) 

RETURN 

END 
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